Spatial Prediction with INLA

Introduction

This vignette illustrates the use of INLA for spatial prediction using examples from Blangiardo and Cameletti (2015) and Illian, Sørbye, and Rue (2012). For prediction of continuous spatial processes, the Lindgren, Rue, and Lindström (2011) stochastic partial differential equations (SPDE) approach is used to approximate the process through an areal Gaussian Markov random field (GMRF) representation. Finally, Log-Gaussian Cox process models are fit using the pseudodata approach of Simpson et al. (2016).

GMRF Background

Blangiardo and Cameletti (2015) section 6.1.

• Observations aggregated to disjoint areal regions indexed by \(i\).
• Each region has unique parameter \(\theta_{i}\).
• \(\mathcal{N}(i)\) is the set of indices of neighbors of region \(i\) and \(\mathcal{N}_{i} = |\mathcal{N}(i)|\) is the number of neighbor of region \(i\).
• Local Markov property: given \(\boldsymbol{\theta}_{\mathcal{N}(i)}\), \(\theta_{i}\) is independent of all other \(\theta_{j}\).

• Then the precision matrix \(\mathbf{Q}\) of \(\boldsymbol{\theta}\) is sparse because only neighbors have nonzero coprecisions.

• Besag-York-Molliè model:
  • Exchangeable random effects \(u_{i}\) with intrinsic conditional autoregressive (iCAR) structure.
  • \(u_{i}|\mathbf{u}_{-i} \sim \mathrm{N}\left(\mu_{i} + \sum_{j} a_{ij} (u_{j} - \mu_{j}) / \mathcal{N}_{i}, \sigma_{u}^{2} / \mathcal{N}_{i}\right)\)
  • (What are the \(a_{ij}?\)).
  • iCAR is an improper prior because covariance matrix not positive definite but this is ok for random effects.

SPDE Background

• Spatial process \(\xi(\mathbf{s})\).

• Solve \((\kappa^{2} - \Delta)^{\alpha / 2}(\tau \xi(\mathbf{s})) = \mathcal{W}(\mathbf{s})\).
  • \(\kappa\) is a range parameter.
  • \(\Delta\) is the Laplacian.
  • \(\alpha\) is a smoothness parameter.
  • \(\tau\) a precision parameter.
  • Exact and sationary solution: \(\xi(\mathbf{s})\) is a Gaussian field with Matèrn covariance function.

• Finite element approximation…

Geostatistics Example

Toy dataset from Blangiardo and Cameletti (2015).

# Plot the data.
plot(s2 ~ s1, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0), data = SPDEtoy, pch = 19, asp = 1, main = 'Toy Data')

# Create a mesh for the SPDE method and then plot it.
toy_mesh <- inla.mesh.2d(as.matrix(SPDEtoy[,c('s1', 's2')]), max.edge = c(0.1, 0.2))
plot(toy_mesh, asp = 1)
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)

# SPDE projector matrix for estimation.
A_est <- inla.spde.make.A(toy_mesh, as.matrix(SPDEtoy[,c('s1', 's2')]))

# Initialize exponential covariance structure for SPDE.
spde <- inla.spde2.matern(mesh = toy_mesh, alpha = 2)

# Set up stack for estimation.
stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = spde$n.spde)
stack_est <- inla.stack(data = list(y = SPDEtoy$y), A = list(A_est), effects = list(c(stack_index, list(intercept = 1))), tag = 'est')

# Create a grid for prediction.
toy_nx <- 50
toy_ny <- 50
toy_grid <- expand.grid(x = seq(0, 1, length.out = toy_nx), y = seq(0, 1, length.out = toy_ny))

# SPDE projector matrix for prediction.
A_pred <- inla.spde.make.A(mesh = toy_mesh, loc = as.matrix(toy_grid))

# Set up stacks for prediction.
stack_latent <- inla.stack(data = list(xi = NA), A = list(A_pred), effects = list(stack_index), tag = 'pred_latent')
stack_response <- inla.stack(data = list(y = NA), A = list(A_pred), effects = list(c(stack_index, list(intercept = 1))), tag = 'pred_response')

# Join all three stacks.
stacks <- inla.stack(stack_est, stack_latent, stack_response)

# Fit the model with INLA.
toy_fit <- inla(
  y ~ -1 + intercept + f(spatial_field, model = spde),
  data = inla.stack.data(stacks),
  control.predictor = list(A = inla.stack.A(stacks), compute = TRUE)
)

# Output posterior summaries.
toy_fit$summary.fixed

toy_fit$summary.hyperpar

# Extract posterior mean of latent spatial field.
index_latent <- inla.stack.index(stacks, tag = 'pred_latent')$data
post_mean <- toy_fit$summary.linear.predictor[index_latent, 'mean']
post_sd <- toy_fit$summary.linear.predictor[index_latent, 'sd']

# Plot the posterior mean and SD of the latent spatial field.
plot(im(matrix(post_mean, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior Mean of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)

plot(im(matrix(post_sd, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior SD of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)

Bei Dataset

Example from Møller and Waagepetersen (2007), Beilschmiedia pendula Lauraceae locations in a plot in Panama. bei dataset in spatstat (Baddeley and Turner 2005).

# Plot the full point pattern.
plot(bei, pch = '.', cols = 'black', main = 'Realized Point Pattern')

# Take a sample of quadrats and plot the observed point pattern.
set.seed(84323)
N_QUADS <- 10
QUAD_SIZE <- 50

w_edge <- Frame(bei)$xrange[1]
e_edge <- Frame(bei)$xrange[2]
s_edge <- Frame(bei)$yrange[1]
n_edge <- Frame(bei)$yrange[2]

botleft <- cbind(
  runif(N_QUADS, w_edge, e_edge - QUAD_SIZE),
  runif(N_QUADS, s_edge, n_edge - QUAD_SIZE)
)
bei_interior <- lapply(seq_len(nrow(botleft)), function(r){return(
    cbind(
      botleft[r, 1] + c(0, 0, QUAD_SIZE, QUAD_SIZE),
      botleft[r, 2] + c(0, QUAD_SIZE, QUAD_SIZE, 0)
    )
  )})
bei_win <- do.call(
  union.owin,
  apply(botleft, 1, function(x){return(
    owin(x[1] + c(0, QUAD_SIZE), x[2] + c(0, QUAD_SIZE))
  )})
)
bei_hole <- bei[complement.owin(bei_win, frame = Frame(bei))]
bei_samp <- bei[bei_win]
bei_window_full <- Window(bei)

plot(bei_hole, main = 'Observed Region with Holes', pch = '.', cols = 'black')

plot(bei_window_full, main = 'Observed Sample')
plot(bei_win, add = TRUE)
plot(bei_samp, pch = '.', cols = 'black', add = TRUE)

# Parameters to experiment with.
MAX_EDGE_LENGTH <- 25
MAX_EDGE_EXT <- 50
MARGIN <- 100

# Mesh covering the site.
bei_boundary <- inla.mesh.segment(loc = do.call(cbind, vertices.owin(Window(bei))))
bei_full_mesh <- inla.mesh.create(
  boundary = bei_boundary,
  refine = list(max.edge = MAX_EDGE_LENGTH)
)
plot(bei_full_mesh, asp = 1)
points(bei, pch = 19, cex = 0.25, col = 'red')

# Mesh including a margin outside the site.
margin_mesh <- inla.mesh.2d(
  loc = bei_full_mesh$loc[,1:2], # Include nodes from site.
  offset = MARGIN,
  max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
margin_spde <- inla.spde2.matern(mesh = margin_mesh)
plot(margin_mesh, asp = 1)
points(bei, pch = 19, cex = 0.25, col = 'red')

# Meshs with coarser resolution in quadrats.
quad_hole <- do.call(
  inla.mesh.segment,
  lapply(seq_along(bei_interior), function(i){
    return(inla.mesh.segment(loc = bei_interior[[i]], grp = i - 1))
  })
)
bei_hole_mesh0 <- inla.mesh.create(
  boundary = list(bei_boundary, quad_hole),
  refine = list(max.edge = MAX_EDGE_LENGTH)
)
plot(bei_hole_mesh0, asp = 1)
points(bei_hole, pch = 19, cex = 0.25, col = 'red')

bei_hole_mesh <- inla.mesh.create(
  loc = bei_hole_mesh0$loc[,1:2], # Include nodes from mesh with holes.
  boundary = bei_boundary,
  refine = list(max.edge = MAX_EDGE_EXT) # Fill in the rest with a coarser triangulation.
)
bei_hole_spde <- inla.spde2.matern(mesh = bei_hole_mesh)
plot(bei_hole_mesh, asp = 1)
points(bei_hole, pch = 19, cex = 0.25, col = 'red')

# Meshs with finer resolution in quadrats.
quad_bnd <- do.call(
  inla.mesh.segment,
  lapply(seq_along(bei_interior), function(i){
    return(inla.mesh.segment(loc = apply(bei_interior[[i]], 2, rev), grp = i - 1))
  })
)
bei_samp_mesh0 <- inla.mesh.create(
  boundary = quad_bnd,
  refine = list(max.edge = MAX_EDGE_LENGTH)
)
plot(bei_samp_mesh0, asp = 1)
points(bei_samp, pch = 19, cex = 0.25, col = 'red')

bei_samp_mesh <- inla.mesh.create(
  loc = bei_samp_mesh0$loc[,1:2], # Include nodes from mesh in quads.
  boundary = bei_boundary,
  refine = list(max.edge = MAX_EDGE_EXT) # Fill in the rest with a coarser triangulation.
)
bei_samp_spde <- inla.spde2.matern(mesh = bei_samp_mesh)
plot(bei_samp_mesh, asp = 1)
points(bei_samp, pch = 19, cex = 0.25, col = 'red')

# Meshes with varying resolution in quadrats and a margin.
margin_hole <- inla.mesh.2d(
  loc = bei_hole_mesh$loc[,1:2], # Include nodes from mesh with holes.
  offset = MARGIN,
  max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
margin_hole_spde <- inla.spde2.matern(mesh = margin_hole)
plot(margin_hole, asp = 1)
points(bei_hole, pch = 19, cex = 0.25, col = 'red')

margin_samp <- inla.mesh.2d(
  loc = bei_samp_mesh$loc[,1:2], # Include nodes from quads.
  offset = MARGIN,
  max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
margin_samp_spde <- inla.spde2.matern(mesh = margin_samp)
plot(margin_samp, asp = 1)
points(bei_samp, pch = 19, cex = 0.25, col = 'red')

obs_full <- rep(0, margin_mesh$n)
obs_full[inla.over_sp_mesh(as(Window(bei), 'SpatialPolygons'), margin_mesh, 'vertex')] <- 1

obs_hole <- rep(0, margin_hole$n)
obs_hole[inla.over_sp_mesh(as(Window(bei_hole), 'SpatialPolygons'), margin_hole, 'vertex')] <- 1

obs_samp <- rep(0, margin_samp$n)
obs_samp[inla.over_sp_mesh(as(Window(bei_samp), 'SpatialPolygons'), margin_samp, 'vertex')] <- 1

Bei Dataset with gridding

NGRID_X <- 40
NGRID_Y <- 20

centers <- gridcenters(
  dilation(bei_window_full, max(NGRID_X, NGRID_Y)),
  NGRID_X, NGRID_Y)
dx <- sum(unique(centers$x)[1:2] * c(-1, 1)) / 2
dy <- sum(unique(centers$y)[1:2] * c(-1, 1)) / 2
bei_df <- data.frame(x = centers$x, y = centers$y,
                     count = NA_integer_, area = NA_real_)

message('gridding full: ', system.time(
for(r in seq_len(nrow(bei_df))){
  bei_df$count[r] <- sum(bei$x >= bei_df$x[r] - dx &
                         bei$x < bei_df$x[r] + dx &
                         bei$y >= bei_df$y[r] - dy &
                         bei$y < bei_df$y[r] + dy)
  bei_df$area[r] <- area(Window(bei)[owin(c(bei_df$x[r] - dx, bei_df$x[r] + dx), c(bei_df$y[r] - dy, bei_df$y[r] + dy))])
}
)['elapsed'], ' seconds')

par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(bei_df$count, nrow = length(unique(bei_df$x)))), unique(bei_df$x), unique(bei_df$y), unitname = 'meters'), ncolcours = range(bei_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(bei_window_full, border = 'white', add = TRUE)
points(bei, pch = '.', col = 'black')

# SPDE projector matrix for estimation.
full_A_est <- inla.spde.make.A(
  margin_mesh,
  as.matrix(bei_df[bei_df$area > 0, c('x', 'y')])
)

# Set up stack for estimation.
stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = margin_spde$n.spde)
stack_est <- inla.stack(data = list(count = bei_df$count[bei_df$area > 0], larea = log(bei_df$area[bei_df$area > 0])), A = list(full_A_est), effects = list(c(stack_index, list(intercept = 1))), tag = 'est')

# SPDE projector matrix for prediction.
full_A_pred <- inla.spde.make.A(mesh = margin_mesh, loc = as.matrix(bei_df[,c('x', 'y')]))

# Set up stacks for prediction.
stack_latent <- inla.stack(data = list(xi = NA), A = list(full_A_pred), effects = list(stack_index), tag = 'pred_latent')
stack_response <- inla.stack(data = list(count = NA), A = list(full_A_pred), effects = list(c(stack_index, list(intercept = 1))), tag = 'pred_response')

# Join all three stacks.
stacks <- inla.stack(stack_est, stack_latent, stack_response)

# Fit the model with INLA.
message('gridded full: ', system.time(
bei_full_fit <- inla(
  count ~ -1 + intercept + f(spatial_field, model = margin_spde),
  offset = larea, family = 'poisson',
  data = inla.stack.data(stacks),
  control.predictor = list(A = inla.stack.A(stacks), compute = TRUE),
  verbose = TRUE
)
)['elapsed'], ' seconds')

# Output posterior summaries.
bei_full_fit$summary.fixed

bei_full_fit$summary.hyperpar

# Extract posterior mean of latent spatial field.
index_pred <- inla.stack.index(stacks, tag = 'pred_latent')$data
post_mean <- bei_full_fit$summary.linear.predictor[index_pred, 'mean']
post_sd <- bei_full_fit$summary.linear.predictor[index_pred, 'sd']

# Plot the posterior mean and SD of the latent spatial field.
plot(im(t(matrix(post_mean, nrow = length(unique(centers$x)), ncol = length(unique(centers$y)))), unique(centers$x), unique(centers$y)), main = 'Posterior Mean of Spatial Field')
plot(bei_window_full, add = TRUE)
points(bei, pch = '.', col = 'black')

plot(im(t(matrix(post_sd, nrow = length(unique(centers$x)), ncol = length(unique(centers$y)))), unique(centers$x), unique(centers$y)), main = 'Posterior SD of Spatial Field')
plot(bei_window_full, add = TRUE)
points(bei, pch = '.', col = 'black')

beihole_df <- data.frame(x = centers$x, y = centers$y,
                         count = NA_integer_, area = NA_real_)

message('gridding holes: ', system.time(
for(r in seq_len(nrow(beihole_df))){
  beihole_df$count[r] <- sum(bei_hole$x >= beihole_df$x[r] - dx &
                             bei_hole$x < beihole_df$x[r] + dx &
                             bei_hole$y >= beihole_df$y[r] - dy &
                             bei_hole$y < beihole_df$y[r] + dy)
  beihole_df$area[r] <- area(Window(bei_hole)[owin(c(beihole_df$x[r] - dx, beihole_df$x[r] + dx), c(beihole_df$y[r] - dy, beihole_df$y[r] + dy))])
}
)['elapsed'], ' seconds')

par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(beihole_df$count, nrow = length(unique(beihole_df$x)))), unique(beihole_df$x), unique(beihole_df$y), unitname = 'meters'), ncolcours = range(beihole_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = '#00000040')

# SPDE projector matrix for estimation.
hole_A_est <- inla.spde.make.A(
  margin_hole,
  as.matrix(bei_df[beihole_df$area > 0, c('x', 'y')])
)

beisamp_df <- data.frame(x = centers$x, y = centers$y,
                         count = NA_integer_, area = NA_real_)

message('gridding sample: ', system.time(
for(r in seq_len(nrow(beisamp_df))){
  beisamp_df$count[r] <- sum(bei_samp$x >= beisamp_df$x[r] - dx &
                             bei_samp$x < beisamp_df$x[r] + dx &
                             bei_samp$y >= beisamp_df$y[r] - dy &
                             bei_samp$y < beisamp_df$y[r] + dy)
  beisamp_df$area[r] <- area(Window(bei_samp)[owin(c(beisamp_df$x[r] - dx, beisamp_df$x[r] + dx), c(beisamp_df$y[r] - dy, beisamp_df$y[r] + dy))])
}
)['elapsed'], ' seconds')

par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(beisamp_df$count, nrow = length(unique(beisamp_df$x)))), unique(beisamp_df$x), unique(beisamp_df$y), unitname = 'meters'), ncolcours = range(beisamp_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = '#00000040')

# SPDE projector matrix for estimation.
samp_A_est <- inla.spde.make.A(
  margin_samp,
  as.matrix(bei_df[beisamp_df$area > 0, c('x', 'y')])
)

Bei Dataset with Simpson et al. (2016) method

This method relies upon the Lindgren, Rue, and Lindström (2011) approximation of the latent Gaussian field as a linear combination of a finite number of basis functions represented as a GMRF on the nodes of a triangulation of the space. Simpson et al. (2016) use the triangulation for numerical integration of the intensity function and show that the LGCP likelihood factors into the joint distribution of independent Poisson random variables corresponding to the points of the point pattern and the nodes of the triangulation. The model fitting proceeds using INLA to fit a Poisson model to pseudodata.

The pseudodata are constructed as follows.

• Let \(n\) be the size of the point pattern.
• Let \(p\) be the number of nodes of triangulation.
• Define the pseudodata as \(\mathbf{y} = (y_{1}, \dots, y_{p}, y_{p+1}, \dots, y_{p+n})'\) where \(y_{i} = 0\) for \(i = 1, \dots, p\) (the traingulation nodes) and \(y_{i} = 1\) for \(i = p+1, \dots, p+n\) (the observed points).
• Define \(\boldsymbol{\alpha} = (\alpha_{i}, \dots, \alpha_{p}, \alpha_{p+1}, \dots, \alpha_{p+n})'\) to encode the numerical integration scheme where \(\alpha_{i}\) is the numerical integration weight for \(i = 1, \dots, p\) (the traingulation nodes) and \(\alpha_{i} = 0\) for \(i = p+1, \dots, p+n\) (the observed points).

Then \(y_{i} \sim Poisson(\alpha_{i}\eta_{i})\) where \(\log(\eta_{i})\) is the SPDE representation of the GF at the location of the \(i\)th pseudodatum. See the paper for tedious notation regarding the definition of \(\eta_{i}\). Ultimately, the nodes become Poisson random variables with means equal to the intensity at that their respective locations, observed points become Poisson random variables with means of 1, and the likelihood is approximately proportional to

$ {i=1}^{p+n} {i}^{y_{i}} (-{i} {i}). $

(Is there a missing \(\alpha_{i}\)?)

NPIX_X <- 400
NPIX_Y <- 200

full_pts <- cbind(bei$x, bei$y)

# Contruct the SPDE A matrix for nodes and points.
full_nV <- margin_mesh$n
full_nData <- dim(full_pts)[1]
full_LocationMatrix <- inla.mesh.project(margin_mesh, full_pts)$A
full_IntegrationMatrix <- sparseMatrix(i = 1:full_nV, j = 1:full_nV, x = rep(1, full_nV))
full_ObservationMatrix <- rbind(full_IntegrationMatrix, full_LocationMatrix)

# Get the integration weights.
full_IntegrationWeights <- diag(inla.mesh.fem(margin_mesh)$c0)
full_E_point_process <- c(obs_full * full_IntegrationWeights, rep(0, full_nData))

# Create the psuedodata.
full_fake_data <- c(rep(0, full_nV), rep(1, full_nData))

# Fit model to full site.
full_formula <- y ~ -1 + intercept + f(idx, model = margin_spde) # No covariates.
full_data <- list(y = full_fake_data, idx = 1:full_nV, intercept = rep(1, full_nV))
message('pseudodata full: ', system.time(
result_full <- inla(
  formula = full_formula,
  data = full_data,
  family = 'poisson',
  control.predictor = list(A = full_ObservationMatrix),
  E = full_E_point_process,
  verbose = TRUE
)
)['elapsed'], ' seconds')
result_full$summary.fixed

result_full$summary.hyperpar

# Plot surface.
proj_margin_mesh <- inla.mesh.projector(margin_mesh, dims = c(NPIX_X, NPIX_Y))
plot(im(t(inla.mesh.project(proj_margin_mesh, result_full$summary.random$idx[,'mean'])),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior Mean Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

plot(im(t(inla.mesh.project(proj_margin_mesh, result_full$summary.random$idx[,'sd'])),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior SD Log-Intensity')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

hole_pts <- cbind(bei_hole$x, bei_hole$y)

# Contruct the SPDE A matrix for nodes and points.
hole_nV <- margin_hole$n
hole_nData <- dim(hole_pts)[1]
hole_LocationMatrix <- inla.mesh.project(margin_hole, hole_pts)$A
hole_IntegrationMatrix <- sparseMatrix(i = 1:hole_nV, j = 1:hole_nV, x = rep(1, hole_nV))
hole_ObservationMatrix <- rbind(hole_IntegrationMatrix, hole_LocationMatrix)

# Get the integration weights.
hole_IntegrationWeights <- diag(inla.mesh.fem(margin_hole)$c0)
hole_E_point_process <- c(obs_hole * hole_IntegrationWeights, rep(0, hole_nData))

# Create the psuedodata.
hole_fake_data <- c(rep(0, hole_nV), rep(1, hole_nData))

# Fit model to site with holes.
hole_formula <- y ~ -1 + intercept + f(idx, model = margin_hole_spde) # No covariates.
hole_data <- list(y = hole_fake_data, idx = 1:hole_nV, intercept = rep(1, hole_nV))
message('pseudodata hole: ', system.time(
result_hole <- inla(
  formula = hole_formula,
  data = hole_data,
  family = 'poisson',
  control.predictor = list(A = hole_ObservationMatrix),
  E = hole_E_point_process,
  verbose = TRUE
)
)['elapsed'], ' seconds')
result_hole$summary.fixed

result_hole$summary.hyperpar

# Plot surface.
proj_margin_hole <- inla.mesh.projector(margin_hole, dims = c(NPIX_X, NPIX_Y))
plot(im(t(inla.mesh.project(proj_margin_hole, result_hole$summary.random$idx[,'mean'])),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior Mean Log-Intensity')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

plot(im(t(inla.mesh.project(proj_margin_hole, result_hole$summary.random$idx[,'sd'])),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior SD Log-Intensity')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

samp_pts <- cbind(bei_samp$x, bei_samp$y)

# Contruct the SPDE A matrix for nodes and points.
samp_nV <- margin_samp$n
samp_nData <- dim(samp_pts)[1]
samp_LocationMatrix <- inla.mesh.project(margin_samp, samp_pts)$A
samp_IntegrationMatrix <- sparseMatrix(i = 1:samp_nV, j = 1:samp_nV, x = rep(1, samp_nV))
samp_ObservationMatrix <- rbind(samp_IntegrationMatrix, samp_LocationMatrix)

# Get the integration weights.
samp_IntegrationWeights <- diag(inla.mesh.fem(margin_samp)$c0)
samp_E_point_process <- c(obs_samp * samp_IntegrationWeights, rep(0, samp_nData))

# Create the psuedodata.
samp_fake_data <- c(rep(0, samp_nV), rep(1, samp_nData))

# Fit model to quadrat-sampled site.
samp_formula <- y ~ -1 + intercept + f(idx, model = margin_samp_spde) # No covariates.
samp_data <- list(y = samp_fake_data, idx = 1:samp_nV, intercept = rep(1, samp_nV))
message('pseudodata sampled: ', system.time(
result_samp <- inla(
  formula = samp_formula,
  data = samp_data,
  family = 'poisson',
  control.predictor = list(A = samp_ObservationMatrix),
  E = samp_E_point_process,
  verbose = TRUE
)
)['elapsed'], ' seconds')
result_samp$summary.fixed

result_samp$summary.hyperpar

# Plot surface.
proj_margin_samp <- inla.mesh.projector(margin_samp, dims = c(NPIX_X, NPIX_Y))
plot(im(t(inla.mesh.project(proj_margin_samp, result_samp$summary.random$idx[,'mean'])),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior Mean Log-Intensity')
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

plot(im(t(inla.mesh.project(proj_margin_samp, result_samp$summary.random$idx[,'sd'])),
        xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
        yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
        unitname = c('meter', 'meters')),
        main = 'Posterior SD Log-Intensity')
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

Bei Dataset and inlabru

bei_full_spdf <- as.SpatialPoints.ppp(bei)
cmp_full <- coordinates ~ mySmooth(map = coordinates, model = margin_spde) + Intercept
message('inlabru full: ', system.time(
bei_full_lgcp <- lgcp(cmp_full, bei_full_spdf, options = list(verbose = TRUE))
)['elapsed'], ' seconds')

bei_full_lgcp$summary.fixed

bei_full_lgcp$summary.hyperpar

# Plot posterior means and posterior sd.
lambda_full <- predict(bei_full_lgcp, pixels(bei_full_mesh), ~ exp(mySmooth + Intercept))
plot(lambda_full)
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

plot(lambda_full['sd'])
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

bei_hole_spdf <- as.SpatialPoints.ppp(bei_hole)
cmp_hole <- coordinates ~ mySmooth(map = coordinates, model = margin_hole_spde) + Intercept
message('inlabru with holes: ', system.time(
bei_hole_lgcp <- lgcp(cmp_hole, bei_hole_spdf, options = list(verbose = TRUE))
)['elapsed'], ' seconds')

bei_hole_lgcp$summary.fixed

bei_hole_lgcp$summary.hyperpar

# Plot posterior means and posterior sd.
lambda_hole <- predict(bei_hole_lgcp, pixels(bei_hole_mesh0), ~ exp(mySmooth + Intercept))
plot(lambda_hole)
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

plot(lambda_hole['sd'])
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

bei_samp_spdf <- as.SpatialPoints.ppp(bei_samp)
cmp_samp <- coordinates ~ mySmooth(map = coordinates, model = margin_samp_spde) + Intercept
message('inlabru quadrats: ', system.time(
bei_samp_lgcp <- lgcp(cmp_samp, bei_samp_spdf, options = list(verbose = TRUE))
)['elapsed'], ' seconds')

# Plot posterior means and posterior sd.
lambda_samp <- predict(bei_samp_lgcp, pixels(bei_samp_mesh0), ~ exp(mySmooth + Intercept))
plot(lambda_samp)
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

plot(lambda_samp['sd'])
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

References

Baddeley, Adrian, and Rolf Turner. 2005. “Spatstat: An R Package for Analyzing Spatial Point Patterns.” Journal of Statistical Software 12 (6): 1–42.

Blangiardo, Marta, and Michela Cameletti. 2015. Spatial and Spatio-Temporal Bayesian Models with R-INLA. Wiley.

Illian, Janine B, Sigrunn H Sørbye, and Håvard Rue. 2012. “A Toolbox for Fitting Complex Spatial Point Process Models Using Integrated Nested Laplace Approximation (Inla).” The Annals of Applied Statistics, 1499–1530.

Lindgren, Finn, Håvard Rue, and Johan Lindström. 2011. “An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73 (4): 423–98.

Møller, J, and RP Waagepetersen. 2007. “Modern Spatial Point Process Modelling and Inference.” Scandinavian Journal of Statistics 34: 643–711.

Simpson, Daniel, Janine B Illian, Finn Lindgren, Sigrunn H Sørbye, and Havard Rue. 2016. “Going Off Grid: Computationally Efficient Inference for Log-Gaussian Cox Processes.” Biometrika 103 (1): 49–70.